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Abstract 

The role of multimode vibrational dynamics in electron transport through single molecule junc- 
tions is investigated. The study is based on a generic model, which describes charge transport 
through a single molecule that is attached to metal leads. To address vibrationally-coupled elec- 
tron transport, we employ a nonequilibrium Green's function approach that extends a method 
recently proposed by Galperin et al. [Phys. Rev. B 73, 045314 (2006)] to multiple vibrational 
modes. The methodology is applied to two systems: a generic model with two vibrational degrees 
of freedom and benzenedibutanethiolate covalently bound to gold electrodes. The results show 
that the coupling to multiple vibrational modes can have a significant effect on the conductance of 
a molecular junction. In particular, we demonstrate the effect of electronically induced coupling 
between different vibrational modes and study nonequilibrium vibrational effects by calculating 
the current-induced excitation of vibrational modes. 

PACS numbers: 85.65.+h, 71.38.-k, 73.23.-b 



I. INTRODUCTION 



Experimental studies of single molecule conductance^ 3 -^ have revealed a wealth of inter- 
esting transport phenomena and have stimulated great interest in the basic mechanisms that 
determine electron transport on the molecular scale JSlSlZl Thereby, effects due to coupling be- 
tween electronic and nuclear degrees of freedom have been of particular interestPThe small 
size and the low mass of molecules may result in strong coupling between electronic and 
nuclear degrees of freedomPEH Vibrational structures in molecular conductance have been 
observed for a variety of different systemsJ 9 * 11 * 12 * 13 * 14 * 15 * 16 * 17 * 1 -^^ Electronic-vibrational cou- 
pling can result in excitation of the vibrational modes of the molecular bridge as well as local 
heating.^ Conformational changes of the geometry of the conducting molecule are possible 
mechanisms for switching behavior and negative differential resistance®'. Furthermore, the 
observation of vibrational structures in conductance measurements allows the unambiguous 
identification of the molecule in the junction. 

These experimental findings have inspired great interest in the theoretical 
modeling and simulation of vibrationally-coupled charge transport in molecular 
junctions™ 2 1 23 1 24 1 25 1 26 1 27 1 28 1 29 1 30 1 31 1 32 1 33 1 34 1 35 ™ A variety of different approaches have been 
used to study the influence of the vibrational degrees of freedom on single molecule conduc- 
tance, including inelastic scattering theory, density matrix approaches and nonequilibrium 
Green's function methods. Green's functions are particularly well-suited to study many-body 
and vibrational nonequilibrium effects. Employing perturbation theory, nonequilibrium 
Green's function methods have been applied in the off-resonant tunneling regime to study, 
e.g., inelastic tunneling spectra)- 25 1 26 1 36 1 In this regime, the effective electronic-vibrational cou- 
pling is typically small and perturbation theory valid. Galperin et al. have recently proposed 
a method that allows the study of vibrational effects in the resonant transport regime.®! This 
method is based on a polaron transformation of the Hamiltonian and employs perturbation 
theory within a self-consistent scheme to solve the equations of motion for the nonequilibrium 
Green's function. 

In the original formulation,®! the method of Galperin et al. is limited to the treatment of 
a single vibrational mode that is coupled to a thermal bath. In this paper, we extend this 
Green's function method to allow the treatment of several vibrational degrees of freedom. 
Moreover, we outline a scheme that allows the calculation of current-induced vibrational 
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excitation and thus to study vibrational nonequilibrium effects. The methodology is applied 
to two different systems: a generic model of a molecular junction with two active vibrational 
modes as well as charge transport through benzenedibutanethiolate covalently bound to gold 
electrodes based on a recently developed^ first-principles model. 



II. THEORY 
A. Model 



To study vibrationally-coupled electron transport through a molecular junction, we con- 
sider a generic tight-binding model,^^ where a single electronic state localized on the 
molecule is coupled to respective states in the left (L) and right (R) leads by tunneling 
matrix elements V4, 

H = eocW e k c W+ Yl ( V k c l c + V k'Jc k ) (2.1) 
fceL,R fceL,R 

+ VL a a) a a a + ^2 KQ a {c ] c - 5) + Uptfpbp + ^2 QaU a pQp. 

a a p a/3 

Here, c\ and are operators that create an electron in the leads and on the molecular 
bridge, respectively, and ej., €q denote the energies of the corresponding electronic states. 



The nuclear degrees of freedom of the molecule are described in Eq. (2.1 ) within the harmonic 
approximation employing the normal modes of the neutral junction in equilibrium at zero 
bias. Thereby, at Q denotes the creation operator for a normal mode of the neutral molecule 
with frequency Q a and Q a = (a a + ajj is the corresponding displacement operator. 

If an external bias is applied to the junction, electrons will be transferred from the 
leads to the junction and vice versa. Accordingly, the potential energy for the nuclei will 
be distorted. The corresponding interaction potential is approximated employing a linear 
expansion in the displacement operator Q a around the equilibrium geometry of the neutral 
junction, resulting in the electronic- vibrational interaction term 'Yj a \aQa{(^c — 5) in Eq. 



(2.1). The charge density, to which this interaction potential is linked, depends on whether 
the molecular state through which the transport takes place is unoccupied in equilibrium 
(5 = 0, corresponding, e.g., to charge transport through the LUMO) or occupied (6 — 1, 
corresponding, e.g., to the HOMO). 40 The corresponding shift in the equilibrium geometry 
of mode a is given by AQ a = 2X a /fl a . 



To account for the effect of relaxation of the vibrational modes, induced by anharmonic 
interactions and coupling to phonons in the electrodes, we adopt a linear response model 
for vibrational relaxation. 333 ^ Within this model, the normal modes of the molecule are 



coupled linearly in Eq. (2.1 ) to a thermal bath of secondary modes. Thereby, bp denotes the 
creation operator for a bath mode with frequency up and U a p determines the system-bath 
coupling strength. All properties of the bath, which influence the dynamics of the system 
are characterized by the spectral density^ 

J a (uj) = ^U^H^-Up). (2.2) 

P 

In the applications considered below, we have used an Ohmic spectral density 

J a (uj) = ri a Lue-" /Wca . (2.3) 

Here, the characteristic frequency u ca defines the maximum of the spectral density and the 
overall coupling strength is determined by r] a . Both values may depend on the specific 
system mode a. 

To apply (self-consistent) perturbation theory within the nonequilibrium Green's func- 
tions approach considered below, it is expedient to remove the direct coupling terms be- 
tween electrons and vibrations in the Hamiltonian, \ a Q a {c*c — 5). To this end, we apply 

a standard canonical transformation, also referred to as Lang-Firsov (or small polaron) 
transformational^ 

H = e s He~ s 1 (2.4) 

= eoc+c + e * c K + J2 ( V k Xc l c + V kX ] c ] Ck) 
fceL,R fceL,R 

+AtW a A + BtW b B + QlUQ b , 

with X = exp (zPjjW^A) and S = — iP) 1 W^ 1 A(c^c — 5). For notational convenience we 
introduce the vectors (Q a ) Q , = Q a , (P a ) Q = — «(Oa — a a)j = a a , (A) a = X a , the matrices 
(W a ) m , = 5 aa 'Q a , (U) a/3 = U a p, as well as respective quantities for the bath modes. It 
is noted that, in the presence of a bath, the decoupling requires that the eigenvalues of 
the matrix (4UW b " 1 U^W a ^ 1 ) are smaller than unityP corresponding to weak system-bath 
coupling. 
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The electronic- vibrational coupling manifests itself in the transformed Hamiltonian (2.4) 
in a polaron shift of the electronic energy 

e - AtW^A, 5 = 0, 
eo - 6 = { (2.5) 
eo + AtW^A, 5 = 1, 

and in the shift generator X that dresses the molecule-lead couplings 



B. Nonequilibrium Green's function approach 

To describe transport properties for the model introduced above, we employ a nonequilib- 
rium Green's function method, which was proposed by Galperin et alP^ Here, we generalize 
this method for applications to multiple vibrational system modes and outline how this 
method can be used to calculate vibrational properties. 

The central quantity in nonequilibrium Green's function theory is the electronic Green's 
function on the molecular bridge 

G(t,t') = -i{T c c{r)c\r')) H , (2.6) 
= -z(T c c(r)ct(r')X(r)Xt(r')) F , 

where T c denotes the time-ordering operator along the Keldysh contour^ and the subscripts 
H/H indicate the Hamiltonian that is used in the calculation of the respective expectation 
value. Most expectation values considered below refer to H, for which the corresponding 
subscript is omitted in the following. 

Following Galperin et al.j^the electronic Green's function G(t,t') is factorized, 

G(r, t>) » G c {r, t') <T c X(r)Xt(r')> , (2.7) 

into a correlation function of the shift generator, (^Y c X{t)X^{t')^j, and an electronic Green's 
function, G c {r, r'), with 

G c (t,t') = -z(T c c(t)J(t')) 7T . (2.8) 

This factorization is valid in the limit of weak molecule-lead coupling corresponding to a 
relatively long residence time of the electron on the molecular bridge. This is the regime, 
where vibrational effects are expected to be particularly pronounced 
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Based on the Hamiltonian H and employing the equation of motion for the electronic 
Green's function, the following equation for G c (t, t') is obtained 331351 

G c (t,t') = G° c (r,r') + Jdndr 2 G° c (r, ri)E c (n, t 2 )G° c (t 2 , t'), (2.9) 

where G° c denotes the electronic Green's function for vanishing electronic coupling (i.e. = 



0). The self energy S c (r, r'), introduced in Eq. (2.9), comprises all interactions of the 
electronic degrees of freedom on the molecule with those in the leads and the vibrational 
modes and is given by 

£ c (r,r') = \V k \ 2 g k (r,r')(T c X(r')X^r)), (2.10) 

fceL,R 

= £°(r,r')<T c X(r')Xt(r)), 

with gt = — 2(T c Cfc(r)cJ.(r')). The above expression holds to second order in V&, i.e. for weak 
molecule-lead coupling. 

To extend the validity to moderate coupling, which is particularly important to describe 
vibrational effects in resonant electron transport, higher order terms are taken into account. 
To this end, a self-consistent scheme is introduced by replacing the last G° c in the integral 
kernel of Eq. Q by G C W M 



G c (T,r') = G C (T,T') + ydT 1 dT 2 G°(T,T 1 )S c (T 1 ,T 2 )G c (T 2 ,T , ). (2.11) 

It is noted that Eq. ( |2.11[ ) gives the exact electronic Green's function for vanishing electronic- 
vibrational coupling (i.e. X a — 0). 



Projection of Eq. (2.11 ) onto the real time axis, according to analytic continuation rules,' 



and Fourier transformation of the resulting equations gives a Dyson equation 

Gl(E) = G°/ + G°/(E)Xl(E)Gl(E), (2.12a) 

and a Keldysh equation 

G<{E) = Gl(E)Z<(E)G a c (E), (2.12b) 

where we have introduced the retarded (G T C ), advanced (G*), and lesser (Gf) Green's func- 
tions and the corresponding self energies. It should be emphasized that the derivation of the 
compact form of these equations requires time-translational invariance and the existence of 
a steady-state transport regime. 



So far, we have only considered the electronic part of the Green's function G(t,t'). 
The solution of the equations requires also the calculation of the shift generator correlation 
function (T c X(t)X^(t')\. Using a cumulant expansion up to second order in the vibronic 
coupling parameters \ Q /£l a yields^^ 

<T c X(r)Xt(r')) = exp f^W" 1 (D(r, r') - D(r, r)) W^A] , (2.13) 

where D(r, r') = — i (T c P a (r)P) 1 (r / )) denotes the correlation matrix of the momentum 
vector P a , in the following referred to as the vibrational Green's function. In the limit 



Vk — > 0, Eq. (2.13 ) constitutes an analytically exact expression. Employing again an equation 
of motion approach, one can derive Dyson/Keldysh equations for the vibrational Green's 
function D(r, r'). The resulting self energy term H e i(r, r') describes interactions between 
electrons and vibrations. This self energy was obtained by Galperin et al. for a single 
vibrational modeP^ Employing the same level of approximations, we obtain for multiple 
vibrational modes 

n el (r, 7-0 = -zW^AAtW; 1 (S c (r, t')G c {t' , r) + £ c (r', r)G c (r, r')) , (2.14) 

and 

D r (£) = TP* + T>°*(E)I% l (E)iy(E), (2.15a) 
D<(£) = D r (£)n<(£)D a (£). (2.15b) 

The off-diagonal elements of the self energy matrix II e i describe interactions between 
different vibrational modes mediated by the electronic degrees of freedom. The matrix D° 
describes momentum correlations of the vibrations in thermal equilibrium and plays the 
same role for the vibrational Green's function D as G° c for the electronic Green's function 
G c . In the presence of a thermal bath, the matrix D° need not to be of diagonal form, 
especially if the bath degrees of freedom couple the various vibrational modes with each 
other. 

The equations for the electronic and vibrational Green's functions outlined above have to 
be solved self-consistently. We employ the following self-consistent schemed As a starting 
point we use the free electronic and vibrational Green's functions. The free vibrational 



Green's function D° enters the shift generator correlation function according to Eq. (2.13). 



Next, the shift generator correlation function is convoluted with the bare self energy E°, 



which gives the dressed self energy S c according to Eq. (2.10). The dressed self energy E c , 
which now contains interactions of the electronic degrees of freedom on the molecule with 
both the leads and the vibrations, is inserted into the electronic Green's function G c , Eqs. 



(2.12), and into the self energy term n e i of the vibrational Green's function, Eq. (2.14). 



The new Green's function D is obtained from the Dyson/Keldysh equations (2.15), and 
enters the shift generator correlation function. Thus, the self-consistent cycle closes, and all 
steps can be repeated with the updated shift generator correlation function. Convergence is 
reached as soon as the variation of the electronic occupation number, n c = lm[Gf(t = 0)], 
between subsequent iteration steps falls below a threshold of 10~ 7 . This way, we obtain 



self-consistent solutions of Eqs. (2.12) and (2.15). 



C. Observables of interest 



Several observables can be considered to study the effects of vibrational motion on charge 
transport through single molecule junctions. Here, we will focus on the current-voltage 
characteristic, the differential conductance and the vibrational nonequilibrium distribution. 

Based on the self-consistent result for the Green's function, the current through lead K 
(K e {L, R}) induced by an external dc-bias $ is obtained employing the formula of Meir, 
Wingreen and JauhcJSEZlMI 

/K = I / ^7 Kk(£)G c >(£) " KAE)Gf(E)} , (2.16) 

with 

X c ,k(t,t') = ^|\4| 2 ^(r,r')<T c X(r')Xt(r)>, (2.17) 
fceK 

= S° K (r,r')<T c X(r')Xt(r)). 

Thereby, the factor two accounts for spin degeneracy. The differential conductance is given 
by g = d//d$. 

It is noted that the scheme outlined above conserves the number of electrons and thus 
obeys Kirchhoff's law, Jl = — Ir- Furthermore, as there is no direct electron-vibrational 
coupling term in the Hamiltonian H, the symmetrized current / = 1/2 (Jl — ^r) can be 
expressed in terms of a transmission function T(E) 

I = y / ^ (ME) - f R (E)) T(E), (2.18) 
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with 



T{E) = r^+r^i) 2 {Gl{E) " G:m ■ (2 ' 19) 

Here, we have introduced the nonequilibrium distribution function 

Im [E< K (£7)1 

/k(^) = p C ; K F \ J , (2-20) 

and the width function 

r K (£) = -2Im[S r CiK (E)]. (2.21) 

Another interesting observable to investigate vibrationally-coupled electron transport in 
molecular junctions is the nonequilibrium vibrational distribution. Due to current-induced 
excitation and deexcitation of the vibrational modes, the vibrational distribution in the 
stationary state may differ significantly from its equilibrium distribution. To study such 
vibrational nonequilibrium effects, we consider in this work the average occupation number 
of different vibrational modes 

n a = {ai<x a ) H . (2.22) 

Within the self-consistent Green's function approach employed here, n a is given (up to 
second order in the system-bath interaction U a p) by the expression 

1\ / 1\ A 2 \ n c 6 = 0, 

//„ ( A a + - ) Im[(D(t = 0))J -[B a + -)+^{ (2.23) 

l-n c , 5 = 1, 



with 

T2 



B « = H (, 2^2^ + 2iV B (^))- (2.24b) 
Thereby, n c = (c f c) H denotes the stationary population of the molecular electronic state. 



The derivation of Eq. (2.23) is outlined in the Appendix. 
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III. RESULTS AND DISCUSSION 



The methodology outlined above has been applied to two different systems: A generic 
model for a molecular junction including two vibrational modes and charge transport through 
benzenedibutanethiolate coupled to gold electrodes. In the latter system the four most 
strongly coupled vibrational modes have been taken into account. 



A. Generic model system with two vibrational degrees of freedom 

First, we consider a generic model system with a single molecular orbital coupled to two 
vibrational degrees of freedom. The energy of the molecular state is chosen as e = leV 
and is located well above the Fermi energy of the leads, ep = OeV. Thus in equilibrium, 
the molecular state is unoccupied (corresponding to 5 = 0). The leads are modelled by a 
one-dimensional tight-binding model with nearest-neighbor coupling constant (3 = 2eV and 
molecule-lead coupling strength v = 0.1 eV. This results in an unstructured semi-elliptic 
conduction band with self energies^ 

Kk(E) = A° (E) - lr° K (E), (3.1a) 
Kk(E) = i&(E)r° K (E), (3.1b) 
Kk(E) = -i(l-tt)(E)r° K (E), (3.1c) 

where the corresponding level-width function reads 

r° K (E) = pm[^(E-^-A^] , (3.2) 

and is the Hilbert transform of T^. Furthermore, denotes the Fermi distribution in 
the leads, 

l + exp(^j 

A t L(R) is the chemical potential in the leads, and $ the bias voltage. We assume the bias 
voltage to drop symmetrically at the right and the left contact, /il(r) = cf ± $/2. In the 
results presented below, we have used a temperature of k-^T = lmeV. This temperature is 
low enough to study vibrational features undistorted by thermal fluctuations {k^T <^ f2i(2)). 
The parameters of the two vibrational modes of the model are given in Tab. [Tj Each mode 
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frequency 


vibronic coupling 


system-bath coupling 


Ol = O.lOeV 
fi 2 = 0.25 eV 


Ai = 0.06 eV 
A 2 = 0.15 eV 


Tfr = 0.001 

T}2 = 0.001 



TABLE I: Vibrational parameters of the model system with two vibrational modes. 



is coupled to an Ohmic bath as described by Eqs. (2.2), (2.3). Thereby, the characteristic 
frequencies of the bath spectral density were chosen to coincide with the frequency of the 
respective system mode, i.e. u ca = £l a - A relatively weak system-bath coupling strength 
was used, r] a = 0.001, corresponding to vibrational relaxation times of about 0.1-lps. In 
principle, the thermal bath couples the two vibrational modes with each other. This inter- 
action is neglected in the calculation presented below. Hence, the retarded projection of the 
correlation matrix D° has a diagonal form 



6aa 'u 2 -m --ion 



a J - L bath,a 



(3.4) 



with 



- 21m [n 



bath, a 



UJ 



2-kJJuj) 



(3.5) 



The approach presented above can be used to describe two different regimes of electron 
transport: inelastic electron tunneling in the off-resonant regime and inelastic resonant elec- 
tron transport. Fig. [TJshows the conductance as a function of bias voltage in the inelastic elec- 
tron tunneling regime, $ < 1 V. In this off-resonant transport regime, electronic- vibrational 
coupling manifests itself in steps of the conductance at voltages $ = fli, 2fii, Q2, •• that cor- 
respond to the opening of inelastic channels. These channels are related to the excitation of 
vibrational quanta and will be referred to as emission channels in the following. The relative 
step heights reflect the respective transition probabilities determined by the corresponding 
Franck-Condon factors. Absorptive channels, corresponding to deexcitation of vibrational 
quanta, play only a minor role in this regime. This is due to the low temperature and the 
small electric current, which is insufficient to drive the vibrational system far from equi- 
librium. The comparison of the results obtained with (black line) and without (gray line) 
vibrational self energy II e i shows that in the off-resonant regime vibrational nonequilibrium 
effects manifest themselves in a shift of the conductance steps. 
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0.1 0.2 0.3 0.4 

bias voltage O [V] 

FIG. 1: Conductance g for the model with two vibrational modes in the inelastic electron tunneling 
regime. The solid gray line refers to a calculation with vibrations in thermal equilibrium, II e i = 0, 
while the solid black line presents results including nonequilibrium effects, II c i ^ 0. The thin 
dashed lines indicate the voltages corresponding to the opening of vibrationally inelastic channels. 

We next consider the resonant transport regime, which is found for voltages $ > 1.5 V. 
Fig. [2] shows the current and the conductance in this regime based on three different cal- 
culations: full vibronic calculations with (solid black lines) and without (solid gray lines) 
vibrational nonequilibrium effects as well as results of a purely electronic calculation (i.e. 
\ a = 0, black dashed lines). The purely electronic calculation exhibits a single peak in the 
conductance at a bias voltage of $ « 2eo when the molecular level enters the conductance 
window between the chemical potentials of the left and the right lead (/xl > eo > A*r)- The 
subsequent decrease of the current results from the finite width of the conduction bandP^l 

The coupling to the vibrational degrees of freedom manifests itself in steps in the current 
and as peaks in the conductance. For the present model (5 = 0), these resonance structures 
can be approximately associated to transitions between the vibrational states of the neutral 
molecule and those of the molecular anion. Thereby, transitions from vibrational states 
of the neutral molecule to those of the anion dominate in the low-voltage regime, where 
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electronic 
vibronic / no n e i 
vibronic / full n e i 



1.25 1.5 



1.75 2. 2.25 
bias voltage $ [V] 



2.5 2.75 




1.75 2. 2.25 
bias voltage [V] 



FIG. 2: Current (top) and conductance (bottom) for the model with two vibrational modes in 
the resonant transport regime. The dashed black line refers to a purely electronic calculation. The 
solid gray line depicts results with vibrations in thermal equilibrium, II e i = 0, while the solid black 
line presents results including nonequilibrium effects, II c i / 0. 
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1.25 1.5 1.75 2. 2.25 2.5 2.75 
bias voltage $ [V] 

FIG. 3: Population of the electronic state, n c = Im[G<(t = 0)], for the model with two vibrational 
modes. 

the electronic state on the molecular bridge is essentially unoccupied (cf. Fig. [3]). For 
higher voltages, the electronic population on the molecular bridge is no longer negligible 
and, therefore, also transitions from the vibrational states in the molecular anion to those in 
the neutral system contribute to the current. Because, in the present model, the vibrational 
frequencies are the same for both electronic states, features related to these transitions 
appear at the same voltages and cannot straightforwardly be distinguished. The population 
of the electronic state is depicted in Fig. [3j The electronic population increases with bias 
voltage in a step-like fashion similar to the current- volt age characteristic (cf. Fig. |2ja)). In 
contrast to the current, which decreases for large voltages due to the finite band width, the 
population increases for larger voltages and exceeds the wide-band limit, which is 0.5 for a 
symmetric junction. 

If vibrational nonequilibrium effects are neglected (II e i = 0, solid gray line), the molecule 
is (due to the low temperature) in its vibrational ground state. Correspondingly, the trans- 
mitting electrons can only induce transitions starting from the vibrational ground state. 
The lowest peak in the conductance corresponds predominantly to the transition from the 
vibrational ground state in the neutral molecule to the vibrational ground state of the anion. 
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Due to the polaron shift (cf. Eq. (2.5)), which accounts for the difference between the vertical 
(purely electronic) and adiabatic transition energy, this peak appears at a lower bias voltage 
($ « 2e , with e = 0.874 eV) than the peak in the purely electronic calculation. Accord- 
ing to the moderate coupling parameters, Ai(2)/fii(2) = 0.6, this conductance peak has the 
largest intensity. Several resonance structures appear at larger voltages. These structures 
can be approximately associated with vibrationally excited states in the molecular anion. 
The relative peak heights follow qualitatively the respective Franck-Condon factors. How- 
ever, the intensities do not coincide with the relative Franck-Condon factors because both 
electron and hole transport contribute to the current. 

If vibrational nonequilibrium effects are included, n e i 7^ 0, additional peaks appear in 
the conductance. These are due to the fact that the nonequilibrium stationary state of 
the vibrational modes is no longer the ground state but involves excited states. Thus, in 
addition to the transitions considered above, where the electron can only lose energy to 
the vibrations, the electrons may induce transitions from vibrationally higher excited states 
to lower vibrational states corresponding to the absorption of vibrational energy by the 
transmitting electrons. As a result, the current (solid black lines) increases already before 
the shifted electronic state enters the conductance window [//l,/^r]. The four conductance 
peaks, which are seen in the inset of Fig. |2](b), are caused predominantly by such absorptive 
processes. Some of the structures involve both excitation and deexcitation of vibrational 
modes. For example, the peak at $ ~ 1.45V is associated to emission of energy to mode (1) 
and absorption of energy from mode (2). Similar processes are found for higher voltages. The 
position of features that exist only due to absorptive processes is highlighted by thin dashed 
lines. For voltages $ > 1.75 V, the current including vibrational nonequilibrium effects is 
found to be smaller than the current neglecting such effects. This is in agreement with 
previous model studies.^ The reduced current reflects the fact that transitions between the 
neutral molecule and the anion are suppressed for a vibrationally excited molecular bridge. 

The findings discussed above are corroborated by Fig. [4], which shows the nonequilibrium 
vibrational excitation 711(2) in the stationary state. For voltages in the resonant regime ($ > 
1.5 V), the results exhibit pronounced vibrational excitation, corresponding to a significant 
deviation of the nonequilibrium vibrational distribution from the equilibrium distribution. 
Because excitation of mode (1) requires less energy than that of mode (2), and because 
of the equal coupling parameters Ai(2)/^i(2) — 0.6, ni is systematically larger than n 2 . 
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1.25 1.5 1.75 2. 2.25 2.5 2.75 
bias voltage O [V] 



FIG. 4: Nonequilibrium vibrational excitation for the model with two vibrational modes. Solid gray 
lines refer to the mode with Oi = O.lOeV, and the solid black lines to the one with O2 = 0.25 eV. 

The vibrational occupation numbers exhibit steplike changes at the same bias voltages as 
the current.^ In contrast to the current- volt age characteristic, however, at some steps the 
vibrational excitation decreases, corresponding to absorptive processes. 

Finally, we discuss the coupling of the two vibrational modes mediated by the electronic 
degrees of freedom. Technically, this coupling is described by the off-diagonal elements of 
the self energy n e i, which depend both on the electronic molecule-lead coupling and the 
electronic-vibrational coupling. Because the coupling of the vibrational modes is mediated 
by the electronic degrees of freedom, it would enter a perturbative description only in second 
order. As a result, for the relatively small molecule-lead coupling considered here, large 
contributions are only expected for quasidegenerate vibrational modes, fli ~ f2 2 . To study 
this effect, Fig. [5] shows the nonequilibrium vibrational excitation of the two modes for fixed 
voltage ($ = 3.5 V) as a function of the frequency of mode (1). Thereby, the electronic- 
vibrational coupling strength Ai was adjusted to keep the dimensionless coupling parameter 
at the constant value X\/fli = 0.6. All other parameters of the model remain unchanged. In 
addition to the results including the full self energy n e i (solid lines), results neglecting the 
off-diagonal elements of II e i are depicted (dashed lines). The results demonstrate that for 

16 



— full n el 

diag(riei) 



0.175 0.2 



0.225 0.25 0.275 
fti [eV] 



0.3 0.325 



FIG. 5: Vibrational excitation numbers versus the vibrational frequency fl± for the model 
with two vibrational modes. The bias voltage is set to 3.5V. Gray lines refer to mode (1), black 
lines to mode (2). Dashed lines are calculated neglecting the off-diagonal elements of n e i, while 
the solid lines are obtained using the full self energy matrix. 



vibrational modes with significantly different frequencies, as considered above in Figs. (2j^4j, 
the electronically mediated mode-mode coupling has little effect and thus the off-diagonal 
elements of II e i could be neglected. For quasidegenerate modes with similar frequencies 
Qi pa Q2, on the other hand, neglecting the off-diagonal elements of II e i would result in 
qualitatively incorrect results. In particular, the results including the full self energy matrix 
predict a resonant behavior of the vibrational excitation for quasidegenerate modes, which 
is missed if the off-diagonal elements of the self energy n e i are neglected. 

It is noted, that if the modes are exactly degenerate, fii = O2 = and have the same 
vibronic coupling, Ai = A 2 = A, the system can be equivalently represented by a single 
vibrational mode with energy Q and coupling parameter a/2A (or yNX for N degenerate 
modes). Calculating the excitation number n singlc for this single mode system, we obtain 
Single ~ ni+n 2 = 2ni(2) . The agreement of n singlc with 2ni(2) corroborates the approximations 



employed in the derivation of Eq. (2.23). 
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frequency 


vibronic coupling in state A 


vibronic coupling in state B 


system-bath coupling 


n a = 0.070 eV 


Ai A) = 0.021 eV 


Ai B) = 0.049 eV 


r} & = 0.001 


O b = 0.149eV 


A b A) = 0.052 eV 


A b B) = 0.037eV 


r/ h = 0.001 


Q c = 0.153eV 


A^ A) = 0.039 eV 


A^ B) = 0.080 eV 


r/ c = 0.001 


Q A = 0.208 eV 


A^ A) = 0.120eV 


aJ, b) = 0.093 eV 


r] d = 0.001 



TABLE II: Vibrational parameters for the benzenedibutanethiolate molecular junction. 



B. Benzenedibutanethiolate 

As a second application, we consider electron transport through p-benzene- 
di(butanethiolate) (BDBT) bound to gold electrodes. This system was chosen because the 
butyl spacer group acts as an insulator between the electronic 7r-system of benzene and the 
gold electrodes. Compared to benzenedithiolatej^ the residence time of the electron on the 
molecular bridge is thus significantly longer, which results in pronounced vibrational effects. 
In recent work, we have developed a first-principles model to describe vibrationally-coupled 
electron transport through BDBT and studied conductance properties employing inelastic 
scattering theory.^ Here, we apply the nonequilibrium Green's function approach outlined 
above to investigate charge transport through BDBT. 

The details of the model are described in Ref. EHJ Briefly, electron transport through 
BDBT is dominated by two electronic states localized at the molecular bridge (denoted A 
and B in the following). The energies of these two states are located at e\ = — 1.38 eV, 
e-Q = — 1.77eV with respect to the Fermi energy of the junction. As a result, the states are 
occupied in equilibrium, corresponding to 5 = 1. The model includes the four most strongly 
coupled vibrational modes. The parameters of these modes are given in Tab. [XTJ Vibrational 
relaxation is described in the same manner as for the model system considered above. 

Strictly speaking, the nonequilibrium approach outlined above can only be applied to a 
single electronic state on the molecular bridge. The extension of the approach to allow the 
description of multiple electronic states will be the subject of future work. In the system 
considered here, the electronic states are well separated from each other (with respect to 
the corresponding level broadening, Tk), and, furthermore, at least three vibrational quanta 
are required to bridge the electronic energy gap, €a — £b ~ 0.37eV. Therefore, we neglect 
coherences between the two electronic states and calculate the overall current as the sum of 
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FIG. 6: Conductance of a benzenedibutanethiolate molecular junction in the inelastic tunneling 



regime. The solid gray line refers to a calculation with vibrations in thermal equilibrium, II e i = 0, 



while the solid black line presents results including nonequilibrium effects, II c i ^ 0. 



the currents obtained separately for the two states. This approximate treatment is supported 
by purely electronic as well as inelastic scattering theory calculations that include both 
electronic states simultaneously.^ 

We start our discussion with the conductance in the inelastic tunneling regime, shown in 
Fig. [6] The results exhibit distinct steps at the onset of the inelastic channels corresponding 
to the excitation of a single vibrational quantum, $ = fi a , f^, ••• For larger voltages, struc- 
tures related to double excitation processes can be seen. These are highlighted with thin 
dashed lines. The comparison of the results with and without self energy II e i shows that the 
renormalization of the vibrational energies, originating from interactions with the electronic 
degrees of freedom, is negligible. 

Next, we consider charge transport through BDBT in the resonant regime. Fig. [7] depicts 
results of calculations with and without vibrational nonequilibrium effects. In addition, 
results of purely electronic calculations (A Q = 0) are shown. The purely electronic result 
(dashed line) shows two pronounced steps at voltages where the two electronic states enter 
the conductance window, respectively. All other smaller structures are due to the energy 
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dependence of the electronic self energy.^ 

The inclusion of electronic-vibrational coupling results in a shift of the two major steps 
from $ = 2.76V to $ = 2.54V and from $ = 3.54V to $ = 3.29V. This shift corresponds 
to the nuclear reorganization energy in the two electronic states. Furthermore, a number 
of additional structures appears. For the model considered here, which is dominated by 
hole transport through occupied states (5 = 1), these structures can be associated with 
transitions from vibrational states in the neutral molecule to those in the molecular cation. 
The first four peaks after the onset of the current can be associated with transitions from 
the vibrational ground state in the neutral molecule to singly excited vibrational states in 
the cation (i.e. 'emission' processes) corresponding to state A. The fourth of these peaks 
shows a significant broadening. We devote this broadening to two mixed emission and 
absorption processes with an energy transfer of AE = Q<i ± (Q c — Q^). For higher voltages, 
a number of structures can be seen, which are related to vibronic transitions in state B. 
Thereby, excitation of single vibrational quanta is the dominant process and gives rise to 
the fundamental peaks at voltages $ = 3.43V, 3.59V, 3.60V, 3.71V. Each of these peaks 
is accompanied by another peak that results from excitation of an additional quantum of 
mode (a) at $ = 3.57V, 3.73V, 3.74V, 3.85V. This result demonstrates the strong coupling 
of mode (a) to state B (cf. Tab. [n]). Absorptive channels, related to transition from higher 
excited vibrational states to lower vibrational states, play only a minor role. Three smaller 
peaks can be assigned to such processes: The first one coincides with the double emission 
peak at $ = 3.57V and corresponds to emission of a vibrational quantum with energy fi^ 
and absorption of a quantum with energy f2 a , i.e. AE = — The second peak at 
$ = 3.46 V involves emission of a vibrational quantum with energy Q c and absorption of 
a quantum with energy Q a , i.e. AE = Q c — f2 a . Finally, the small feature at $ = 3.87V 
corresponds to a transition with AE = + Q c — f2 a . 

Electronically mediated mode-mode coupling described by the off-diagonal elements of 
the self energy matrix n e i has negligible influence on the results (data not shown). In view 
of the very similar frequencies of modes (b) and (c), this finding is at first glance surprising. 
The difference of their frequencies, Q c — fib, is, however, large compared to the electronic 
level broadening T K thus effectively reducing the mode-mode coupling. 

In a recent study, we have investigated vibrational effects in charge transport through 
BDBT employing inelastic scattering theory.^ Inelastic scattering has the advantage that 
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FIG. 7: Current and conductance of the benzenedibutanethiolate molecular junction in the reso- 
nant transport regime. The dashed black lines represent results of a purely electronic calculation 
(A Q = 0). Solid gray lines refer to a calculation with vibrations in thermal equilibrium, II c i = 0, 
while solid black lines present results including nonequilibrium effects, II e i / 0. 
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FIG. 8: Comparison of different methods for the current through a benzenedibutanethiolate molec- 
ular junction. Shown are results of vibronic calculations employing the nonequilibrium Green's 
function approach with the full self energy matrix II c i (solid line) and inelastic scattering theory 
(dotted line) The dashed line depicts results of a purely electronic calculation, where both meth- 
ods give identical results. In contrast to Fig. [7j all calculations include only one electronic state 
(state B). 

the transmission process of a single electron is described numerically exactly. The calcula- 
tion of the current, however, involves the assumption that the vibrational degrees of freedom 
of the molecular bridge relax to their equilibrium distribution between two consecutive elec- 
tron transmission processes. Furthermore, the scattering theory approach cannot account 
straightforwardly for the nonstationary electronic occupation of the molecular states and, 
therefore, fails if the molecular levels, which determine the transport, are located close to 
the Fermi energyP^SD pjg [3] shows a comparison of results obtained for the current through 
BDBT employing scattering theory and the nonequilibrium Green's function approach. In 
addition, results of a purely electronic calculation are shown, where both methods give iden- 
tical results. To allow a direct comparison, the results have been obtained including only 
one electronic state on the molecule (state B) and without coupling to the vibrational bath 
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3.5 4. 4.5 5. 5.5 

bias voltage [V] 

FIG. 9: Nonequilibrium vibrational excitation in charge transport through benzenedibutanethio- 
late. Shown is the average number of vibrational quanta in the stationary state for modes (a)-(d). 
The calculation includes only transport through state B and was obtained for a system-bath cou- 
pling strength of rf a = 1CT 3 . 

(rja = 0). Overall, the results show a rather good agreement between the two methods. This 
is due to the fact that the molecule-lead coupling is rather weak and the energy of state B is 
located well below the Fermi energy. The major differences are an earlier onset of the current 
and a smaller current for larger voltages in the Green's function calculation, as well as dif- 
ferent heights of the step structures. The earlier onset of the current as well as the different 
step heights are caused by absorptive vibrational transitions related to the nonequilibrium 
vibrational distribution and contributions from electron transport. Such processes are not 
accounted for in scattering theory, which assumes the vibrational degrees of freedom to relax 
to their equilibrium distribution between two consecutive electron transmission events and 
(in this case) only considers hole transport. Because the transmission probability is nor- 
malized to unity, both results approach each other for larger biases ($ ~ 3.7V). Increasing 
the bias voltage even further, the results start to deviate again. In this regime, current- 
induced vibrational excitations cannot be neglected and lead to a substantial suppression of 
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FIG. 10: Nonequilibrium vibrational excitation in charge transport through benzenedibutanethi- 
olate. Shown is the average number of vibrational quanta in the stationary state for mode (a) 
for different system-bath coupling strengths, ry a . The calculation includes only transport through 
state B. 

the current in the Green's function calculation.^ 

To conclude this section, we consider nonequilibrium vibrational excitation induced by 
the current through the BDBT junction. Fig. [9] shows the average vibrational excitation 
of the four modes. It is seen that in the resonant transport regime the current results in 
significant deviations from the equilibrium distribution. Overall, the amount of excitation of 
the four modes follows the value of the dimensionless vibronic coupling X a /Q a . In particular 
mode (a), which has the strongest vibronic coupling, acquires pronounced excitations in the 
stationary state. 

Current-induced excitation of vibrational modes competes with relaxation of vibrational 
energy due to system-bath coupling and thus depends on the timescale of vibrational relax- 



ation. The influence of the relaxation time is illustrated in Fig. 10 for the strongest coupled 
mode (a). The results obtained for different values of the system-bath coupling strength, 
r] a , show the increase of vibrational excitation for longer vibrational relaxation times. In the 
limit of negligible vibrational relaxation (j] a = 0), mode (a) is excited to rather high quan- 
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turn numbers (n a > 20). The corresponding energy is still smaller than typical dissociation 
energies of a C-C or a C-H bond (-Ediss > 3eV). In this regime, however, the harmonic 
approximation ceases to be valid and anharmonicities should be taken into account. 

IV. CONCLUSIONS 

In this paper we have studied the effect of multimode vibrational dynamics on charge 
transport through single molecule junctions. To this end, we have extended a nonequilib- 
rium Green's function method, developed by Galperin et al.J^to treat multiple vibrational 
modes in transport calculations. This method is based on a polaron transformation of the 
Hamiltonian and employs perturbation theory within a self-consistent scheme to solve the 
equations of motion for the nonequilibrium Green's function. In addition to the simulation 
of electronic properties such as the electronic current and the conductance, we have also 
outlined a scheme to calculate the average vibrational excitation in the stationary state 
using the nonequilibrium Green's function approach. 

The methodology has been applied to two examples: a generic model of a molec- 
ular junction with two active vibrational modes as well as charge transport through 
benzenedibutanethiolate covalently bound to gold electrodes based on a first-principles 
model. In both cases, the results show that the electronic-vibrational coupling may have 
significant effects on the current through the molecular junction. The coupling to the vi- 
brational degrees of freedom manifests itself in pronounced structures in the current-voltage 
characteristic and the differential conductance. Moreover, the current-induced excitation of 
vibrational modes may result in a significant deviation of the nonequilibrium vibrational dis- 
tribution from the equilibrium distribution. For modes with similar frequency, mode-mode 
coupling mediated by the electronic degrees of freedom is of importance. The results show 
that in this case the full self energy matrix needs to be taken into account in the calculation. 

The study in this paper was based on a model, which describes the vibrational degrees of 
freedom in the harmonic approximation. This approximation is well suited for small ampli- 
tude motion and low excitation energies. To investigate the influence of higher vibrational 
excitation and the possible dissociation of the molecular bridge, the approach has to be 
extended to include anharmonic effects. This will be the subject of future work. 
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APPENDIX A: AVERAGE NUMBER OF VIBRATIONAL EXCITATIONS 



In this Appendix, we outline the method used to calculate the average vibrational exci- 
tation in the stationary state 

n a = (ala a ) H . (Al) 

In contrast to the population of the electronic state n c , the excitation number n a of mode 
a is not an invariant of the Lang-Firsov transformation, 

n c = ( c t c ) H = <ctXXt c )_ = ( c t c )_ = Im[G<(t = 0)] , (A2) 

., , . \-4 a a)n - (QaC j c) 1T + ^tn c , 5=0, 
a a a ) H - \ ^ 2 . . . ^ A6 J 



n. 



*)h - k (Qa^c) ir + ^(l-n c ) + ^ {Q a )„, 8=1. 



As a result, we have to calculate four different expectation values to obtain the average 
excitation number n a . Thereby, the subscripts H/H denote the Hamilton-operator that 
is used to compute the respective expectation values. In the following we consider only 
expectation values calculated with H, for which the corresponding subscript is omitted from 
now on. 



i A 2 

The third term in Eqs. (A3), oc , which can be interpreted as the contribution from 



polaron formation, can directly be extracted from the electronic Green's function G c . It 
only contributes in the resonant transport regime, where n c deviates significantly from its 
equilibrium value S. 



The fourth term in Eqs. (A3) for the case 5=1 contains a single displacement operator 
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Q a . For its determination we consider the steady state relations 

= = -*W a (Q a > - 2zU (Q b ) , (A4a) 

= = -*W b (Q b ) - 2zUt (Q a ) , (A4b) 

which result in 

(1 - 4W a - 1 UW b 1 U+) (Q a ) = 0. (A5) 
Because the eigenvalues of the matrix 4W a " 1 UW b 1 have to be smaller than unity in 



order to perform the Lang-Firsov transformation (cf. Sec. [IT]), the expectation value of Q c 
vanishes, (Q a ) = 0. 



To treat the second term in Eq. (A3), we consider similar steady-state relations 



= ^ 1 = -iW a <Q a ct c ) - 2zU <Q b ct c > 

~ E Vk { V - C l CX ) + H V k <PaC f C fc Xt) ; ( A6a ) 

feeL.R ' fceL,R 

= \ h f 1 = -iW b <Q b ct c ) - 2zUt <Q aC t C ) 

- J2 v * ( p * c i cX ) + E v * (Pb^c^t) . (A6b) 

Here, additional terms linear in the molecule-lead coupling Vk appear. However, because all 
expectation values are taken at equal times, these terms cancel each other. This finding is 
closely related to the existence of a steady-state transport regime that preserves Kirchhoff 's 
law, Jl = —Ir- Without the terms linear in Vk, we can follow the same line of argument as 
for (Q a ), and obtain \Q a c^c) = 0. 

Finally, to calculate ^a a a a ^, we exploit the momentum correlation matrix D: 

(at a a ) = \ {(Q a Q a ) - Im[(D(i = 0)) J) - \ (A7) 

The expectation value (Q a Q a ) can be rewritten in terms of (P a P a ) an d {QpQp)- For that 
purpose, we use the following steady-state relations 

= - {Q f a) = in a (P a P a ) - itt a (Q a Q a ) -2iJ2 U*p (Q a Qp) , (A8a) 
= - {Q d f^ = in a (P a P p ) - iup (Q a Q p ) -2iJ2 U a >p {Q a Q a ') , (A8b) 

a 1 

= td{P g Q t Q ^ = ttup (P a P p ) - tn a (Q a Qp) -2iJ2 U af 3> (QpQp) , (A8c) 

13' 
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where, based on the same argument as before, we have disregarded terms linear in Using 



these relations, we obtain to second order in the system-bath coupling (0(U%g)) 



(QaQa) - (PaPa) + £ O r, ^- C)2 \ { P <* P <«) ~ £ , ,2 ~ZcH (QpQp) 



(A9) 



(pp)(i + V 



/9 W "a 



Thereby, we disregard off-diagonal contributions (P a P a i), which turn out to be negligible. 
This step is corroborated by considering the special case of degenerate modes (cf. the dis- 



cussion at the end of Sec. Ill A). 



Inserting Eq. (A9) into Eq. (A7), we obtain the final expression for the excitation number 



n r 



1\ / 1\ I 2 |n c , 5 = 0, 

A « + 2 J M( D (* = °)) a J ~ [ Ba+ 2J + TE . 

1 n ci 5 1 , 



(A10) 



with 



u 2 



p^p 



u 2 



(All) 



The bath modes remain in thermal equilibrium, for which (QpQp) = 1 + 2Nb(up) with 
Nb(uj) the Bose distribution function. 
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